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Abstract 

A self-consistent field method is developed, which can be used to construct models of differentially 
rotating stars to first post-Newtonian order. The rotation law is specified by the specific angular 
momentum distribution j(m U7 ), where m TO is the baryonic mass fraction inside the surface of 
constant specific angular momentum. The method is then used to compute models of the nascent 
neutron stars resulting from the accretion induced collapse of white dwarfs. The result shows that 
the ratios of kinetic energy to gravitational binding energy, (5, of the relativistic models are slightly 
smaller than the corresponding values of the Newtonian models. 
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I. INTRODUCTION 



In a recent paper |IJ (hereafter Paper II), we demonstrate that the accretion induced 
collapse (AIC) of a rapidly rotating white dwarf can result in a rapidly rotating neutron 
star that is dynamically unstable to the bar-mode instability, which is the instability re- 
sulting from non-axisymmetric perturbations with angular dependence e ip . Here (p is the 
azimuthal angle. This instability could emit a substantial amount of gravitational radiation 
that could be detectable by gravitational wave interferometers, such as LIGO, VIRGO, GEO 
and TAMA. 

However, for this instability to occur, the neutron star must have a (3 = T/\W\ greater 
than a critical value (3d ~ 0.25 (Paper II). Here T is the rotational kinetic energy and \W\ is 
the gravitational binding energy. Only the AIC of white dwarfs that are composed of oxygen, 
neon and magnesium (O-Ne-Mg white dwarfs) with Q > 0.93Q m can produce neutron stars 
with such a high value of (3. Here Q is the angular velocity of the white dwarf and Q m is 
the angular velocity at which mass shedding occurs on the equatorial surface. This type of 
source will not be promising for LIGO II because its event rate is not expected to be very 
high. 

Neutron stars are compact objects. General relativistic effects have a significant influence 
on both the structure and dynamical stability of the stars. Recently, Shibata, Baumgarte, 
Saijo and Shapiro studied the dynamical stability of differentially rotating polytropes in 
full general relativity || and in the post-Newtonian approximation ||]. They performed 
numerical simulations on the differentially rotating polytropes with some specified rotation 
law. They found that as the star becomes more compact, the critical value (3d slightly 
decreases from the Newtonian value 0.26 to 0.24 for their chosen rotation law. It is not 
clear, however, whether their result implies that relativistic effects would destabilize the 
stars we are studying, for the equilibrium structure of the star will also be changed by 
relativistic effects. The value of /3 of a relativistic star will not be the same as that of a 
Newtonian star with the same baryon mass and total angular momentum. 

The objective of this paper is twofold. First, we develop a new numerical technique 
to construct the equilibrium structure of a rotating star with a specified specific angular 
momentum distribution to first post-Newtonian (1PN) order [i.e., including terms of order 
c -2 higher than the Newtonian terms in the equations of motion]. Then we use this new 
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technique to construct models of neutron stars corresponding to the collapse of the white 
dwarfs we studied in Paper II and Ref. |2j (hereafter Paper I) and compare them with the 
Newtonian models. 

Equilibrium models of neutron stars in full general relativity have been built by many 
authors 0, || [j], |8], |], |Tt], |TTJ. The neutron stars studied in the literature are either rigidly 
rotating or rotating with an ad hoc rotation law. New-born neutron stars resulting from core 
collapse of massive stars or accretion induced collapse of massive white dwarfs are differen- 
tially rotating Jl], [|, 12, 13, It seems plausible that the rotation laws of these neutron 
stars could be approximated by the specific angular momentum distribution j(m ro ) of the 
pre-collapse stars (see Paper I and Sec. II A). Here m ro is the baryonic mass fraction inside 
the surface of constant specific angular momentum. Equilibrium models of Newtonian stars 



with a specified j(m ro ) have been constructed by many authors [15, 16, 17, 18|. However, 
none of these studies, to our knowledge, has been generalized to include the relativistic 
effects. 

If a rotating axisymmetric star is described by a barotropic equation of state, i.e., the total 
energy density e is a function of pressure only, then there is a constraint on the rotation 
law (see Section II A). This rotational constraint is usually written in the form u°u v = 
F(Q) || W \, where F is an arbitrary function. Here Q is the angular velocity of the fluid 
with respect to an inertial observer at infinity; u° is the time component of the four- velocity 
and u v = u^ip 11 , where (p a is the axial Killing vector field of the spacetime. In the Newtonian 
limit, this constraint reduces to the well-known result that Q is constant in the direction 
parallel to the rotation axis. The major obstacle in the construction of differentially rotating 
relativistic stars is that it is not clear what function F should be used to produce the desired 
specific angular momentum distribution j(m ro ). In the next section, we will reformulate the 
rotational constraint in a way that can be used to impose the rotation law j(m ro ), at least 
in the 1PN calculations. 

The structure of this paper is as follows. In Section II, we give a brief review on the 
full relativistic treatment of rotating relativistic stars and then reformulate the rotational 
constraint imposed by the barotropic equation of state. Next, we use the standard 1PN 
metric and show that the rotational constraint can be solved analytically. We then derive 
the equations of motion determining the structure of a star to 1PN order. In Section [TTI|, 



we generalize the self-consistent field method of Smith and Centrella |23| so that it can be 
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used to compute the structure of a star to 1PN order. In Section [TV], we apply the numerical 
method to construct neutron star models resulting from the collapse of the O-Ne-Mg white 
dwarfs we studied in Paper II and compare them with the corresponding Newtonian models. 
Finally, we summarize our conclusions in Section M. 



II. FORMALISM 



In this section, we first give a brief review on the full relativistic treatment of rotating 
relativistic stars and then reformulate the rotational constraint imposed by the barotropic 
equation of state (EOS) in Section II A. Then we derive the equations of motion determining 



the equilibrium structure of a rotating star to 1PN order in Section [TIB| . Throughout 
this chapter, we use the convention that Greek indices run from to 4, being the time 
component; whereas Latin indices run from 1 to 3 only. A sum over repeated indices is 
assumed unless stated otherwise. The signature of the metric is ( — h ++). 



A. Relativistic hydrodynamics 

We want to construct the nascent neutron stars resulting from the AIC of rotating white 
dwarfs. As in Paper I, we make the following assumptions on the AIC and the collapsed 
stars. 

First, we assume the collapse is axisymmetric. Hence the spacetime, albeit dynamical, 
has an axial Killing vector field (p a . Second, we neglect viscosity and assume a perfect fluid 
stress-energy tensor 

= (e + P)uV + Pg^ , (1) 

where e is the energy density in the fluid's rest frame; P is pressure and w M is the fluid's 
four- velocity, normalized so that u a u a = —1. Third, we assume that the collapsed ob- 
jects can be described by a barotropic EOS, i.e. e = e(P). Fourth, we assume there is no 
meridional circulation in the equilibrium state of the collapsed objects, i.e., the spacetime 
is nonconvective or circular ||30|| . In other words, the fluid's four- velocity can be written as 

/ = u° (t a + ^<p a ) , (2) 

where t a is the timelike Killing vector field of the spacetime of the collapsed star, c is the 
speed of light, u° is the time component of the four- velocity, and Q is the rotational angular 
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velocity measured by an inertial observer at infinity. Finally, we assume that no material is 

ejected from the star during and after the collapse. Hence the total baryon rest mass Mq 

and the total angular momentum J are conserved. 

Let n denote the baryon number density in the fluid's rest frame. It follows from 

the baryon number conservation law V u (nu u ) = and conservation of energy-momentum 

V„T^ = that (see, e.g., Chapter 22 of Ref. |2J) 

de _e + P 
dn 



n 



(3) 



Given a barotropic EOS, the above equation can be integrated, giving 

de' 



me 



n(e ) exp 



(4) 



',o e' + P{e> 

We define the baryonic rest mass density p = nffiB- Here is the average baryon mass, 
defined so that 



lim 



1 



pc z 



(5) 



It follows |19| from the conservation of baryon number V v {nu v ) = 0, conservation of energy- 
momentum V V T^ U = 0, and the existence of an axial Killing vector (p a that 



= | = 0, 



(6) 



where r is the proper time along the fluid particle's worldline and the specific angular 
momentum 



3 



e + P 
pc 



-u 



tp ■ 



(7) 



Here u v = u a f a . In the Newtonian limit, j = Qw 2 , which is the Newtonian expression of 
the specific angular momentum along the rotation axis. Here w is the distance from the 
rotation axis. Following the arguments in Paper I, we conclude that the final neutron star 
should have the same baryon mass Mo, total angular momentum J, and specific angular 
momentum distribution j(m ro ) as the pre-collapse white dwarf. 

In the stationary and axisymmetric spacetime of a relativistic star, the Euler equation 
takes the form |19| 



' ~ —u^Vpu a = V a (ln-u°) — vPu J^ 7 °^ 



(8) 



e + P 

Since the EOS is barotropic, the left side of Eq. @ is a total differential. This imposes a 
constraint on the rotation law: the integrability condition for Eq. (El) is that the rotation 



law must have the form vPu^ = F(Q) [§, 19|, where F is an arbitrary function. In the 
Newtonian limit, this rotational constraint means that Q is constant in the direction parallel 
to the rotation axis. The constraint written in this form is not convenient for our purpose, 
as our rotation laws are specified by the function j(m m ). Hence, we formulate the constraint 
in another way: the integrability condition is that vPu^V a Q is an exact differential. In the 
language of differential forms, we require that u°u v dQ be an exact form. This implies that 
its exterior derivative vanishes: 

d(u°u^dtt) = , (9) 
where d denotes the exterior derivative. 



B. Post-Newtonian approximation 



Following Chandrasekhar [20|1, we split the energy density into two terms: 



n 

e = P c 1 + 3 



(10) 



We adopt the 1PN metric (in Cartesian coordinates) developed by Chandrasekhar, and 



Blanchet, D amour and S chafer [EG, 21, |22| 
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The metric components satisfy the gauge condition 



J2 d i9ij + \ d i (#00 = 0(c 



In this metric, the components of the four-velocity are 













1 -\ 
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(11) 
(12) 
(13) 



(14) 
(15) 



(16) 
(17) 
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where v l = cu l /u° and v 2 = 5ijV t v J . The potentials U and A{ satisfy the elliptic equations 



DjD j U + ^-^-U = AnGp 



If 9 3P N 

1 + - (U + 2v 2 + — 



cf- 



P 



DjEPAi = l6irGp5 ik v k 



(18) 
(19) 



where Dj denotes the covariant derivative compatible with the three dimensional flat-space 
metric, and G is the gravitational constant. We introduce cylindrical coordinates (w, (p, z) 
with d/d(p being the axial Killing vector, and w = \/x 2 + y' 2 . In this coordinate system, 
the velocity vector potential Ai has only one component: A v = xA y — yA x . Let Q = A v /w. 
Then Q satisfies the equation 



To 1PN order, we have 



D j DjQ 



cu°u 



Q 



lQnGpVtw . 



K 



-Hi ( 1 + — 



where 



K = v 2 - AU + — . 

v 



Since vPu^V a VL is an exact differential, we can write 



K 



V a f = wHl 1 + 



(20) 

(21) 
(22) 

(23) 



where / is a scalar function. The rotational constraint (|9]) gives only one nontrivial equation 
for O in a stationary and axisymmetric spacetime. To 1PN order, Eq. (|9|) can be solved 
analytically, giving 



Q(vj, z) = Slo(zcr) + W ^f l ° [K(w, z) - K Q (w)\ 



(24) 

where VLq(vj) = Q(w,0) and K (w) = K(w,0). Eq. (p3|) can then be integrated and we 
obtain, up to an arbitrary additive constant, 



-y^o( ro 

1 



dw' w'Qq^zd') 



where 



fli(zu, z) 



[K{w,z)-K Q {w)\ 



dw' nl(zu') 



zu 



zu'K (zu') + —d^K (w') 



(25) 

(26) 
(27) 
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It is convenient to define an auxiliary function 

2 f P dp ' 

h = c 2 — — . 28 

Jo e{P') + P' V ; 

This quantity is defined only inside the star. The boundary of the star is given by the surface 
h = 0. In the Newtonian limit, h reduces to the specific enthalpy The Euler equation (^j), 
to 1PN order, can be written in integral form: 



t"OJ 

h(zu,z) = I dw' w'Vti(w') — U + C 
Jo 

+ \ ( ^ Q - 2zu 2 Q 2 U + Qv- + h) , (29) 



where C is a constant and all quantities outside the integral are evaluated at (zu, z). 

The structure of the star is determined once a rotation law is given. The rotation law is 
specified by the specific angular momentum distribution function j(m ro ), which is determined 
by the pre-collapse white dwarf (see Paper I). Straightforward calculations using Eqs. (|7j), 
(0), (O), (0), (0), (0), (0), © and © give 



vo 2 Vt, 



1 v 2 P Oi Q 

cM 2 p il v 



(30) 



To compute m ro , the baryonic mass fraction inside the surface of constant j, we first need 
to determine the surfaces on which j is constant. In the Newtonian case, the surfaces of 
constant j are cylinders. This is not true in general in the relativistic case (at least not in 
the coordinate system we are using). Let [zu + r](zu, z)/c 2 , z] denote the surface of constant 
j that intersects the equatorial plane at cylindrical coordinate radius zu. Hence we have 

j[zu + r](zu,z)/c 2 ,z] = j(w,0) , (31) 
r/(ro,0) = 0. (32) 

Expanding the left side of Eq. ( |3~T| ) to 0(c~ 2 ), we obtain 

/ \ 2 j(^^)-joM 

r){zu, z) = -c , , 33 

<9 ro j (w) 

where jo(zu) = j(zu, 0). Using Eq. (|30|) , we obtain 

_ zun £(U -3U + P/p) + zuQ 1 + £{Q) 
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where £(q) = q(zu,z) — q(zu,0). The baryon mass M ro inside the volume V m bounded by 
the surface of constant j is given by 

M m = [ pu^n^dV 
dz 



2tt 



zup [W, z 



/ dzu' zu' p*(zu' , z') 
Jo 



(35) 
(36) 



where is the unit vector orthogonal to the surface of constant t; dV is the proper volume 
element in the constant t hypersurface, and 



P = P 



1 + 



3U 



The baryonic mass fraction is then given by 



M n 



(37) 



(38) 



where Mq is simply the value of M m at zu = R e , the equatorial radius of the star. It is 
convenient to define the normalized specific angular momentum 



Jr 



M . 



J 



Straightforward calculations from Eq. (|30|) give 



■J 



n 2 



ZU H 



a z Q I w 2 n 2 - qu + 2n + 



2-Pq \ 2^oQo 



(39) 



(40) 



Po J ^ 

where A = J/M and the subscript "0" in the above equation means that the quantity is 
evaluated at (zu,0). The integrated Euler equation (E^) becomes 



Here 



h = \ 2 il) - U + C + \ (- — 2Uv 2 + Qv- lv 2 K + I 2 
r \ 4 2 



— dzu , 

o zu'- 5 



dzu' 



v 2 Q d w K Q + 2zu'n 2 U 



Tin 



Po, 



Qo^c 



(41) 

(42) 
(43) 



where all the quantities in the integrands are evaluated at w — w'. 

The rotational kinetic energy T and gravitation potential energy W of a relativistic star 
are given by (see, e.g., |TT|) 

T = \j ndJ =Y c J VT^n^ dV , (44) 
W = -[(M p - M)c 2 + T] , (45) 



where the proper mass M p and gravitational mass M are 



M p = 1 / ,-,/"//„ ,/V f 10) 



M = ~J (T af3 - l -T°g a ^ t a n p dV . (47) 

Both T and W are independent of gauge in a spacetime that is stationary, axisymmetric 
and nonconvective. The expressions for T and W to 1PN order are 



i + \ (v 2 - su + n + - + Q 

C z \ p V 



W = J(pv 2 + pU + 3P) d 3 x-^J p (Jf/ 2 + 



d 3 x , (48) 

Hrr 2 9 4 

— Uv v 

2 8 

Qv „ 2 ^ rr 3Pw 2 6PU\ ,o 

-- — n^ 2 -mi + — d 3 x , 49 

2 2p p J y } 

where d?x = w dw dz dip. 

Given the total baryon mass M , total angular momentum J, normalized specific angular 
momentum distribution j n (m m ), and EOS, we have to solve Eqs. (p~8|) , (p0|), (|4"ID , (p0|). 
and (|38j) consistently to determine the structure of the differentially rotating star. We shall 
discuss how these equations can be solved numerically in the next section. 

III. NUMERICAL METHOD 

In this section, we develop a self-consistent field technique to calculate the structure of a 
relativistic star with the rotation law specified by the normalized specific angular momen- 
tum distribution j n (m ro ). Our method is a generalization of the one used by Smith and 



Centrella [Z4|, which is a modified version of Hachisu's self-consistent field method p5 | . 

The self-consistent field method is an iteration procedure. Suppose in a certain iteration 
step, we have h(w,z) and £lo(m) in a cylindrical grid, we first evaluate the quantities p, 
P and II from the EOS. Then we compute the potentials U and Q by solving the elliptic 



equations ( |i8| ) and fl20[) . Since the velocity potential Q always appears in the 1PN terms of 
the equations of motion, we can replace Q on the right side of Eq. ( p0|) by Qq. The angular 
velocity Q, as well as v = zufl, outside the equatorial plane are determined by Eqs. ( |22|) 
and (^4j). Next, we compute the baryonic mass fraction m m using Eqs. (0), (|36|) , (|37D 
and fl38|) . The function if> is then calculated by Eq. (0). During each iteration, we fix two 
parameters, which we choose to be the central energy density e c [or equivalently, h c = h(0, 0)] 
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and the equatorial radius R e . The constants C and A 2 in Eq. ([|l]) are then given by 

C = h c + U c , (50) 



^ 2 - -7 



U-C-\(i-2v>U + Qv- l ^ + I: 



2 



(51) 



c 2 \A 

where U c = U(0, 0) and all the quantities in the second equation are evaluated at the 



equatorial surface of the star. Finally, we update h by Eq. ( flip and update Qq by solving 
the algebraic equation pp|). We repeat the procedure until h and Q converge to the desired 
accuracy. 

When the star becomes flattened, the iteration scheme described above does not converge. 



This is fixed by generalizing the modified scheme suggested in Ref. ||17|| : the variables h and 
fl in the {i + l)-th iteration, h i+1 and (f2 )i+i are changed to 

h i+1 = hi6 + h'(l - 5) , (52) 
(Q ) i+1 = (n ) t 5 + %(l-6) , (53) 

where h' and Q' are the quantities determined by Eqs. ([41]) an d (0)- ^ ne parameter 5 (0 < 
5 < 1) is used to control the changes of h and Qq in an iteration step. For a very flattened 
configuration, we need to use 5 > 0.9 to ensure convergence, and it takes more than 100 
iterations for the models to converge to a fractional accuracy of 10~ 5 . In the standard self- 
consistent field method, one only needs to solve for the density distribution p (or equivalently, 
the enthalpy distribution h) self-consistently. In our self-consistent field method, we also 
need to solve for the equatorial angular velocity distribution Qq self-consistently. This is the 
main difference between the standard scheme and our proposed scheme, apart from the fact 
that the equations of motion in the 1PN case are more complicated. 

The self-consistent field method described above computes stars with a given central 
energy density e c and equatorial radius R e . However, we want to construct a star with a 
given total baryon mass M and total angular momentum J. To do this, we first compute a 
model of non-rotating spherical star by solving the 1PN structure equations for nonrotating 
stars in isotropic coordinates. We use the density distribution as an initial guess to construct 
a model with slightly different e c and R e . We then build models with different values of e c and 
R e until we end up with the model having the desired baryon mass and angular momentum. 

For a rapidly rotating configuration, the equatorial radius extends to R e > 1000 km and 
the polar radius is approximately R p ~ 10 km in our coordinate system. Hence we use a 
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nonuniform cylindrical grid to perform most of the computations. The resolution near the 
center of the star is about 0.4 km, whereas the resolution is about 6.5 km near the equatorial 
surface of the star. We double the resolution to check the convergence. For a given e c and 
R e , the fractional differences of the baryon mass M and angular momentum J between the 
two resolution grids are less than 10~ 5 even for the rapidly rotating cases. 

We adopt the Bethe- Johnson EOS [EH] for densities above 10 14 g cm" 3 , and BBP EOS 



for densities in the range 10 11 — 10 14 g cm 3 . The EOS for densities below 10 11 g cm 3 is 
joined by that of the pre-collapse white dwarfs, which is the EOS of a zero-temperature 



ideal degenerate electron gas with electrostatic corrections derived by Salpeter ||28|| . We are 
mainly interested in the structure of the most rapidly rotating neutron stars. The central 
densities of these stars are around 4 x 10 14 g cm -3 (see the next section), and ideas about 
the EOS in this relatively low density region have not changed very much since 1970's. 

The baryon masses of the neutron stars we compute in this chapter are around 1.4M . 
For a non-rotating spherical star of this baryon mass, c 2 R/GM « 6 for the EOS we adopt. 
Here M « 1.3M Q is the gravitational mass and R ~ 12 km is the circumferential radius of 
the star. Hence we expect that the second and higher order post-Newtonian terms will give 
about 1/6 2 ~3% corrections to our models. 

The effect of higher post-Newtonian terms can also be estimated by the virial identity to 
1PN order (see the Appendix). We define a dimensionless quantity 

^w N > < 84 > 



where 



Z = J [av 2 + 3P - a(zod m U + zd z U)] d 3 x 

+— J [QPU + pQ zuQ + (2P - av 2 - ApU)(wd m U + zd z U) 
+pwVLo(wd. w Q + zd z Q)} d 3 x , (55) 

W N = - J a{vjd^U + zd z U) d 3 x , (56) 



a 



1 + ^ K - 2U + n + ^ 



(57) 



The dimensionless quantity ( is a measure of the fractional error of our equilibrium models 
due to higher post-Newtonian corrections (see the Appendix). 
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FIG. 1: The central densities e c /c 2 of differentially rotating neutron stars as a function of Q,/Q m 
of the pre-collapse white dwarfs. Both Newtonian and 1PN results are shown for stars having the 
same Mq and J. 

IV. RESULTS 

We only construct neutron star models corresponding to the collapse of O-Ne-Mg white 
dwarfs (i.e., the Sequence III white dwarfs in Paper II), because these neutron stars are the 
most likely to undergo a dynamical instability and emit strong gravitational waves. 

Figure [I] shows the central densities e c /c 2 of the resulting neutron stars as a function 
of Q/Q m , where Q is the angular frequency of the pre-collapse white dwarf, and Q m is the 
angular frequency of the maximally rotating white dwarf in the sequence. Both Newtonian 
and 1PN results are shown for stars having the same Mq and J. We see that the central 
energy densities for the 1PN models are higher than the Newtonian models. This is ex- 
pected because relativistic effects tend to make the stars more compact. The difference in 
e c decreases as the star becomes more rapidly rotating. 

Figure || shows the value of (3 — Tj\W\ of the neutron stars as a function of fl/fl m for 
both the Newtonian and 1PN models. To estimate the possible error of (3 due to higher 
post-Newtonian effects, we use the formula 



We found that the quantity £/ 2 /c 4 is the largest second post-Newtonian terms we neglected 




(58) 
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FIG. 2: The value of ft of the resulting neutron stars as a function of £l/VL m of the pre-collapse 
white dwarfs. The vertical bars through the 1PN curve show /3± 6ft of selected equilibrium models, 
where 5(3 is the error of (3 due to higher post-Newtonian corrections estimated by Eq. (|58|). 



in the whole calculation. Hence we estimate that 

5T 
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5W 



1 

T 



U 2 
Ti—d 3 x 
c 4 



W 



1 r U 2 

w 



(59) 
< 60 ' 

where T, and Wi are the integrands in Eqs. (^) and (|49|), respectively. The vertical bars 
in Fig. |] show ft ± 5/5 of selected equilibrium models [using Eqs. (|58|)-(|60D for 5ft}. The 
result suggests that relativistic effect lowers the value of ft for stars of given M and J. 
The maximum ft of these neutron stars is 0.24, which is 8% lower than the Newtonian case 
(0.26). However, the error bars also suggest that higher post-Newtonian corrections could 
change the values of ft significantly. 

Figure [3] shows the virial quantity ( for our equilibrium neutron star models. This 
quantity is a measure of the fractional correction to the equilibrium structure of the stars 
due to higher post-Newtonian effects. We see that ( is smaller than 0.03 for all the models 



we computed, which is in accord with our estimate near the end of Sec. |Tl I 



The structure of the neutron stars is not much different from the Newtonian models. 
Stars with ft > 0.1 all contain a high-density central core of size about 20 km, surrounded 
by a low-density torus-like envelope. The size of the envelope ranges from 100 km (for stars 
with ft ~ 0.1) to over 500 km (for ft > 0.2). Figure ^ shows the density contour of a typical 
rapidly rotating neutron star. This figure looks basically the same as Fig. 3 of Paper II, 
which shows the density contours of the same star computed with Newtonian gravity. The 
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FIG. 3: The virial quantity £ defined in Eq. (|54j ) for the equilibrium neutron star models, 
parametrized by Q/Q. m °f the pre-collapse white dwarfs. 
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FIG. 4: Meridional density contours of the neutron star resulting from the AIC of a rigidly rotating 
O-Ne-Mg white dwarf with f2/fi m = 0.964. This neutron star has = 0.238. The contours in the 
upper graph denote, from inward to outward, e/e c = 10 -4 , 10 -5 , 10~ 6 , 10~ 7 , 10 -8 , 10 -9 and 0. 
The contours in the lower graph denote, from inward to outward, e/e c = 0.8, 0.6, 0.4, 0.2, 0.1, 10~ 2 , 
10~ 3 and 10 -4 . The central density of the star is e c /c 2 = 3.8 x 10 14 g cm -3 . 

f3 of this star is 0.238, which is somewhat smaller than the Newtonian value 0.255. 

Figure |5] shows the equatorial angular velocity distribution Qq(zj) of the same star. We 
see that the angular velocity in the inner core of the star {w ^$ 20 km) in the 1PN model is 
slightly larger than that of the corresponding Newtonian model. This is expected because 
relativistic effects make the star more compact. The material is compressed more in the 1PN 
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FIG. 5: The equatorial angular velocity Qq(vj) of the neutron star in Fig. |j. 
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FIG. 6: The distribution of the rotational kinetic energy T and gravitational binding energy \W\ 
of the material inside the radius w, for the neutron star in Fig. |j. 

model, and should rotate faster due to the conservation of angular momentum. Figure |6| 
shows the distribution of rotational kinetic energy T and gravitational binding energy \W\ 
of the material contained within cylindrical radius w. The two quantities approach their 
asymptotic values at w ~ 30 km. This is due to the high central condensation of the star. 
Both T and \W\ in the IPN model are larger than the corresponding Newtonian model. The 
kinetic energy T is larger because the star rotates faster. However, the difference between 
the two T-curves decreases as we move away from the rotation axis. This is because most 
of the kinetic energy of the star is from the region 10 km < vj < 30 km, in which relativistic 
effects are less important. On the other hand, the gravitational binding energy is mainly 
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FIG. 7: The equatorial rotational angular velocity fio as a function of vj for w < 60 km. These 
neutron stars result from the AIC of the pre-collapse white dwarfs with Q/Cl m equal to (a) 0.090, 
(b) 0.23, (c) 0.30, (d) 0.41, (e) 0.55, (f) 0.71, (g) 0.86 and (h) 1.00. 




03 (km) 

FIG. 8: The cylindrical mass fraction TYh'pj ELS cl function of w for the neutron star models in Fig. [?]. 
The curves and the corresponding models are identified by the degrees of central condensation: the 
higher the degree of central condensation, the lower the value of 0/O m . 

contributed from the material in the inner region w < 20 km, in which relativistic effects 
are important. As a result, the T/|W| value of the relativistic model is somewhat less than 
in the corresponding Newtonian model. 

Figure [7] shows the equatorial angular velocity Qq for several selected models in the 
central region near the rotation axis. The shape of the curves are very similar to those of 
the Newtonian models (see Fig. 4 of Paper II). However, the angular velocities in the 1PN 
models are all slightly larger than the Newtonian models in the inner core. 



17 




0.25 



g 



f 



c 



d 







10 



100 



1000 



G5 (km) 



FIG. 9: The value of /3 ro as a function of w for the neutron star models in Fig. |7[ 



Figure |^ shows the baryonic mass fraction m TO verses w for the selected models in Fig. [7]. 
As in the Newtonian case (see Fig. 5 of Paper II), the mass is highly concentrated in the 
inner core of the star. The degree of central condensation decreases as the star rotates faster. 
However, more than 80% of the mass is contained within a 30 km radius even for the most 
rapidly rotating star, where the outer envelope extends to over 1000 km. The collapsed 
object can be regarded as a neutron star of size about 20 km surrounded by an accretion 
torus. 

Figure |9| shows the T/|VK| of the material inside the surface of constant w, for the 
selected models in Fig. |7|. The shape of the curves are qualitatively the same as those in the 
Newtonian models ( see Fig. 6 of Paper II), although the values of are slightly smaller. 
All the curves level off at w ~ 30 — 50 km, suggesting that the material in the outer layers 
does not have much influence on the overall dynamical stability of the star. 

V. CONCLUSIONS 

We have generalized the self-consistent field method so that it can be used to compute 
models of differentially rotating stars to 1PN order with a specified angular momentum 
distribution j(m m ). We also applied this new method to construct models of nascent neutron 
stars resulting from the collapse of massive O-Ne-Mg white dwarfs we studied in Paper II 
and compare them with the corresponding Newtonian models. 

We found that the 1PN models are more compact and rotate faster. Our calculations also 
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suggest that the 1PN models have smaller values of (3 than the corresponding Newtonian 
models. The highest value of (3 these neutron stars can achieve is 0.24, which is 8% smaller 
than the Newtonian case. Higher post-Newtonian corrections may change the values of /3, 
but our estimate suggests that they will still be smaller than those of the Newtonian models. 
We estimate that the fractional error of our 1PN models due to our neglecting higher order 
post-Newtonian terms is about 3%. 

The fact that relativistic effects lower the values of (3 can be understood by the following 
heuristic argument. Relativistic effects make a star more compact. Hence the gravitational 
binding enetgy \W\ increases. The compactness of the star also make it rotate faster due to 
conservation of angular momentum. Hence the rotational kinetic energy T also increases. 
However, most of the kinetic energy comes from material slightly away from the center of the 
star (in the region 10 km < w ^$ 30 km for the star in Fig. ||), where relativistic effects are 
less important compared to the region near the center. On the other hand, the gravitational 
binding energy |W| is contributed mainly from material near the center of the star (in the 
region w ^$ 20 km for the star in Fig. |6|). As a result, the increase in T is not as much as 
the increase in and so the value of j3 decreases. 

We have demonstrated that relativistic effects could lower the value of (3 of a star with 
a given baryon mass Mq and angular momentum J. Shibata, Baumgarte, Shapiro and 
Saijo H |J demonstrated that relativistic effects also lower the critical value (3d for the 
dynamical instability by a similar amount. It will be interesting to find out which of these 
two effects is more important. Careful numerical 1PN stability analyses must be carried out 
to determine whether or not relativistic effects destabilize the stars. 
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APPENDIX A: VIRIAL THEOREM TO 1PN ORDER 



Virial theorems in full general relativity were formulated by Bonazzola and Gourgoul- 



hon [29, 31, 32], and have been proved to be useful as a consistency check for numerical 
computation of rotating star models (see, e.g., || |33, 



I. These virial identities are not very 



convenient to implement since they involve integrals over the entire two or three-dimensional 
volume, but our computational domain only extends to the region slightly outside the sur- 
face of the star. Chandrasekhar derived a virial identify to the 1PN order [^U] in which 
integrations are only performed over the interior of a star. However, it involves a double 
volume integral. In this appendix, we shall derive another virial identify which is correct 
to the 1PN order and is easier to compute than the Chandrasekhar identity. This expres- 
sion is useful to estimate the error of our equilibrium neutron star models due to higher 
post-Newtonian effects. 



We start with Eq. (67) of Ref. p0|, which is the Euler equation to 1PN order |35| 



di{av^) + (1 + f j d 3 P + a(l + ^) d 3 U - ^fxvjv'diU 



a 



2f/ + n + 



(Al) 
(A2) 



Here we have set all partial time derivatives to zero. The equation is valid to 1PN order 



for the metric (|TT|)-(|T3|). Multiplying Eq. ( |A1| ) by x 3 , summing over j and integrating over 
d 3 x = dx dy dz = w dw dz dip, we obtain, after integration by parts, 



J {av 2 + 3P - aaPdjU + ^[2P{SU + x j djU) - av 2 x j djU + ^pxhjV%V 



-pv'x^d.Aj - djAi) - ApUxWjU}} d 6 x = 



(A3) 



In a stationary, axisymmetric and circular spacetime, x J Vj = x 3 Aj = 0. In cylindrical 
coordinates, we have 



x 3 djU = tud^U + zd z ll , 
v l x 3 (djAi — diAj) = m£l(zud^Q + zd z Q) + wVlQ 



(A4) 
(A5) 



where Q = (xA y — yA x )/w. Hence Eq. ( |A3| ) becomes, to 1PN order, 
J {av 2 + 3P - a{wdJJ + zd z U) + ^{QPU + P Q wQ 
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+ (2P - av 2 - ApU)(wd m U + zd z U) + pwtt^wd^Q + zd z Q)]}d 3 x = . (A6) 

Equation ( |A6| ) is our virial identity. It involves a volume integral over the interior of the 
star. In the Newtonian limit, it reduces to the familiar form 2T + 3A + W = 0, where 
A = JPd 3 x. 

Let Z be the left side of Eq. ( |A6|) and define a dimensionless quantity 

C = Z/W N , (A7) 
W N = - J o{wdJJ + zd z U)d 3 x . (A8) 

In the Newtonian limit, ( = (2T + 3A + W) /W, which is a quantity often used as a measure 
of the fractional numerical error of Newtonian equilibrium models caused by finite grid size. 

Our equilibrium models are constructed by solving Eq. (|4l|) , which agrees with Eq. (|A1|) 
to 1PN order. Hence the quantity ( measures the fractional error of our models due to 
second and higher post-Newtonian corrections, and also the numerical error due to finite 
grid size. 
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